Universality in the run-up of shock waves 1 
We investigate the run-up of a shock wave from inside to the surface of a perfect 

fluid star in equilibrium and bounded by vacuum. Near the surface we approximate the 

fluid motion as plane-symmetric and the gravitational field as constant. We consider the 

"hot" equation of state P = (V — l)pe and its "cold" (fixed entropy, barotropic) form 

P = Kop T (the latter does not allow for shock heating). We find numerically that the 

evolution of generic initial data approaches universal similarity solutions sufficiently near 

the surface, and we construct these similarity solutions explicitly. The two equations of 

state show very different behaviour, because shock heating becomes the dominant effect 

when it is allowed. In the barotropic case, the fluid velocity behind the shock approaches 

a constant value, while the density behind the shock approaches a power law in space, as 

the shock approaches the surface. In the hot case with shock heating, the density jumps 

by a constant factor through the shock, while the sound speed and fluid velocity behind 

the shock diverge in a whiplash effect. We tabulate the similarity exponents as a function 

of the equation of state parameter T and the stratification index n*. 
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1. Introduction 

In the numerical simulation of high-energy astrophysical events, one is often faced 
with three linked unknowns: the mathematical model of the relevant physics (such as 
an equation of state), the qualitative features to expect (such as shocks) and a numer- 
ical method that correctly represents these features. This is the case for the surfaces of 
compact binaries just before merger. 

The full simulation of compact binary mergers in general relativity has been achieved 
by a small number of groups in recent years. Recent papers describing their methods 



include Baiotti et al. (2008) and Kiuchi et al. (2009) for binary neutron star mergers, and 



Etienne et aL(2009) and Kyutoku et al. (2010) for black hole- neutron star mergers, all 



using high-resolution shock-capturing (HRSC) methods (see Font(2008) for a review) 



for simulating the matter, and Oechslin et al. (2002) an cjBauswein et al. (2010) for binary 
neutron star mergers using smoothed particle hydrodynamics (SPH) methods. 

Given the complications involved in general relativity and the large computational re- 



Universality in the run-up of shock waves 3 
quirements in 3-dimensional simulations involving a range of length scales, in all these 

simulations the matter is modelled as a single perfect fluid. The equations of state used 

range from the two simple families that we also consider here, over composite polytropcs 

and tabulated realistic cold equations of state to combinations of a tabulated cold com- 



ponent and a simple hot component. See Duez(2009) for a review. If we take current 
single fluid models seriously, the neutron star or stars in these models have a surface at 
finite radius where both the density and pressure reach zero, with vacuum outside. (By 
contrast, if a star has a solid crust, the density is small but finite at the surface, and the 
results of this paper do not apply.) 

Both HRSC and SPH methods have serious difficulties in correctly modelling the 
physics in this region. This is essentially because quantities such as the density and pres- 
sure approach zero as powers of distance to the surface, and hence are not smooth at the 
surface itself. This affects in particular the numerical modelling of the balance between 
the gravitational field and pressure gradient in a star in equilibrium. HRSC methods, 
such as those cited above, replace the vacuum exterior by an unphysical low-density 
"numerical atmosphere" with a numerically modified time evolution. SPH methods put 
few particles there. As a consequence, current binary evolution codes are incapable of 
correctly representing behaviour involving shocks so near the surface that the physical 
density there becomes comparable to the atmosphere density. A typical value for the 
density of the numerical atmosphere is 10~ 9 of the central neutron star density, low in 
relative terms, but dense enough in absolute terms, 10 9 kg/m 3 , for a perfect fluid approx- 
imation to hold. 

If one is mainly interested in the generation of gravitational waves by the bulk motion 
of mass, the low-density, low-pressure exterior regions can probably be safely neglected. 
However, photon and neutrino emission are sensitive to temperature as well as mass, 
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and therefore could be dominated by shock heating at the surface. Sound waves running 
down the density and pressure gradient towards the surface have a tendency to shock, 



see 



Gundlach & Please(2010) for a semi-quantitative criterion. In particular, tidal forces 



in a binary will raise waves in the interior of each star, which can form shocks as they 



approach the surface, see Gundlach & Murphy(2010) Very close to the surface, these 
shocks approach similarity solutions. This means that they become singular, but at the 
same time that their behaviour is universal, depending on the initial data only through 
a single parameter. 

As a step towards a correct numerical simulation of such phenomena in astrophysi- 
cal contexts, we construct the similarity solutions as the solution of ordinary differential 
equations, and compute the related similarity exponents. We then use the Clawpack soft- 



ware of LeVeque et aL(2010) to show numerically that the similarity solution is indeed 
an attractor as the shock approaches the surface. To study the self-similarity it is nec- 
essary to zoom in by many orders of magnitude on the origin (the point where vacuum 
is reached) as the computation proceeds. The Clawpack software has been modified to 
repeatedly cut the domain in half and redistribute the computational points over the new 
domain in order to continue the calculation well into the asymptotic regime. Similar ap- 
proaches have previously been used to study self-similarity numerically, e.g. in the study 



of self-similar blow up in Berger & Kohn(1988) Our approach is described in Sec.[6j 

Throughout this paper, we assume that the fluid motion is plane-parallel, and that 
the gravitational field is constant. The former is likely to be a good approximation for 
a sufficiently small part of the surface, sufficiently close to the surface, but this assump- 
tion will need to be tested in two or three-dimensional simulations in future work. The 
approximation of constant gravitational field as g — GM/R 2 is safe when combined 
with the assumption of plane-parallel symmetry because the gravity of the outer layers is 
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neglible compared to the bulk mass M of the star, and the distance R to the centre varies 

little. As we are only considering a small region near the surface of star, general rela- 
tivistic gravity can safely be approximated by Newtonian gravity, and the Lanc-Emdcn 
and Tolman-Oppenhcimer-Volkov solutions for spherical stars agree near the surface. 

We assume perfect fluid matter in the sense of local thermal equilibrium, and consider 
two related equations of state. Throughout this paper, "hot EOS" will refer to the two- 
parameter equation of state P{p, e) = (T — l)pe often called the Gamma-law EOS, and 
"cold EOS" to the one-parameter (barotropic) equation of state P{p) = Kop r often 
called the polytropic EOS. The former reduces to the latter if the entropy is constant. 
The latter does not allow for shock heating (and when we use it we do not solve the 
energy conservation law), and we consider it here mainly for comparison with the more 



physical hot EOS. In Baiotti et aL(2008) these two equations, for T = 2, are compared 
explicitly in neutron star binary mergers. 

In this paper, we consider only Newtonian fluid motion. We find that this is self- 
consistent in the cold case. In the hot case, the fluid velocity and sound speed calculated in 
Newtonian fluid dynamics diverge as the shock approaches the surface, and so, depending 
on the initial data and on how close to the surface the perfect fluid approximation breaks 
down, one may have to go to special relativity. 

Sees. [2] and [3] review the fluid equations and the Rankinc-Hugoniot conditions at the 
shock. 



Seed] generalises the similarity solution of Barker & Whitham(1980) from the shallow 
water case T = 2 to a cold polytropic gas with arbitrary T. We add a discussion of the 
evolution after the shock reaches the surface. 



Sec. [5]rcderives the similarity solution of Sakurai(1960) for the non-relativistic motion 
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with the hot EOS, and tabulates more cases. We add a discussion of the equilibrium 
reached long after the shock has reached the surface. 

In Sec. [6] we present details of the numerical code, and numerical evidence that generic 
smooth perturbations of the hydrostatic equilibrium are attracted to the similarity solu- 
tion as the surface is approached. Sec. [7] presents our conclusions. 

In the Appendix, general mathematical facts about power-law similarity solutions in 
one space and one time coordinate are derived, and then applied fluid motion with our hot 
and cold EOS. We discuss under what circumstances a conservation law in the evolution 
equations gives rise to an integral of the motion in the ODE system obtained under 
a similarity ansatz, and we show how the similarity solutions with the cold equation of 
state relate to the subclass of isentropic similarity solutions with the related hot equation 
of state. 

In the following, ~ denotes equality up to sub-leading terms, while ~ denotes equality 
up to sign, constant factors and sub-leading terms. In particular, our convention is such 
that t and x are usually negative. To simplify notation, we will take care of such signs 
when using the symbol ~, but will neglect them when using ~. The symbol = indicates 
definitions. 

2. The equations 



The ID mass conservation and momentum conservation laws in a constant gravitational 
field g in the —x direction are 



2.1. Conservation laws 



Pt + (vp) x = 0, 



(2.1) 



[vp) t + (v 2 p + P) 



X — 



gp- 



(2.2) 
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Here p is the mass density, v the velocity, and P the pressure. In the cold case, these 



equations arc closed with the barotropic equation of state 



P( P ) = K oP T , 



(2.3) 



Here Kq and T = 1 + 1/n are constants with range Kq > and n > (and hence T > 1). 



The n = 1 cold case is equivalent to the shallow water equations (see LeVeque(2002) I 
and the source term of (|2.2p in this case corresponds to a linear beach with slope 1, so 
that the similarity solution in this case shows the behavior of a bore on a linear beach. 
By contrast, in the hot case the equation of state is 



P(p,e) = (T-l)pe, 



(2.4) 



where the evolution of the internal energy per rest mass e is given by the total energy 
conservation law 



p ( e + — ] v + Pv 



-gpv. 



(2.5) 



2.2. Thermodynamic equations 
The first law of thermodynamics in the form 

de = Tds - Pdip- 1 ), 



(2.6) 



where s is the entropy per rest mass and T the temperature, can be integrated over a 
curve of constant entropy in thcrmodynamical state space, with P given by (|2.4|) to give 



e{P,s)^— -p , 



(2.7) 



where K(s) is some, still unspecified, function of the entropy. Hence we can write the 
hot equation of state in the form 



P(p,s) = K(s)p I 



(2.8) 
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The cold equation of state can therefore be obtained as the isentropic case K = Kq 
const of the hot equation of state with the same polytropic index T (or n). 
This form of the equation of state can also be used to show that 



2 = dP 
dp 



TP 

- (2.9) 



for both the hot and cold case. (Physically, c > is the sound speed.) In the hot case, 
this can also be written as 

c 2 = -e. (2.10) 

n 

Making some additional independent assumption, such as the ideal gas law P = RTp, 
would determine the function K(s), but for the remainder of this paper, we do not need 
to make such a choice. 

2.3. Riemann invariants 
In the cold case, the equations for smooth solutions can be rewritten in the simpler form 

[d t + {v ± c)d x ] (v ± 2nc) + g = 0. (2.11) 

For g = 0, the equations therefore have a pair of Riemann invariants v ± 2nc, with 
characteristic speeds v ± c respectively. In the case g ^ 0, the Riemann invariants still 
exist, and are v + gt ± 2nc, with the same speeds. 

It is interesting to write the hot equations in almost Riemann form as 

[d t + vd x ] K = 0, (2.12) 
[d t + {v± c)d x ] (v ± 2nc) +g = -^c 2 d x In K, (2.13) 

where 

K = Pp- r (2.14) 
is in physical terms a function of only the entropy per mass, and is therefore conserved 
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along particle trajectories. Note that in the isentropic case the hot evolution equations 

reduce to the cold ones. 

2.4. Equilibrium solution 
Under our assumptions, the equation of hydrostatic equilibrium is 

P x = -gp. (2.15) 

In the cold case, the density in hydrostatic equilibrium, after adjusting the origin of x, 
is given by 

P(«) = (~z£rY (2-16) 



nTK, 



which implies 



~f (2.17) 



for x < 0, where x = is the surface of the star. 

In the hot EOS case in general, the hydrostatic equilibrium solution depends on the 
entropy stratification of the star. In the astrophysics literature, a type of stratification is 
often considered where 

P (x) = K*p{x) v % (2.18) 

where and T* are constants which describe the equilibrium star. We define r» = 
1 + l/n». Note that the equation of state is still given by (|2.8p . with K not in general 
constant. Using the hydrostatic equilibrium condition ()2.15|) and the expression (|2.9j) for 
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the sound speed, we find the generalised equilibrium solution 



<x) = _gxn^ (22Q) 

nA * 



/-|^, (2.21) 



1- 



tf(z)=lU ^— ) ■ (2.22) 

(The last two expressions are redundant). The case n* = n reduces to the cold result 
above (with K(x) = A'*). The range n < < oo gives a stratification that is stabilised 
by buoyancy forces against (non-planar) vertical displacements. We will therefore only 
consider n» > n. The marginally stable limit n* = n corresponds to constant entropy, 
something that can be achieved by convection, while = oo can be interpreted as the 
isothermal limit. We will not consider the latter, as the star then has no surface at finite 
radius. 

3. Shock into fluid at rest 

3.1. Mass and momentum conservation 

Let v s , p s , P s be the quantities just behind a shock travelling in the positive x-direction 
(towards the surface at x = 0) with speed er, with v = and po and Pq just in front of 
the shock. Defining 

M - "A (3.1) 
with the range < fx < 1, the mass conservation law gives 

- = 1- M . (3.2) 
(j 

Using (|2.9p . and defining the dimensionless quantity 
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with range < A < oo, the momentum conservation law gives 



1-fi \Po 



11 



(3.4) 



3.2. Cold equation of state 



In the cold case 



and so 



Pn 



ir r , 



A 



(i- At )(i-M r r 

In the limit A — > we have [i — > 0, with 

i 

A s 



/' 



(3.5) 



(3.6) 



(3.7) 



In the opposite limit A — > oo, which will not be relevant here, we have ~ 1 — A 1 ' 2 . 

3.3. Energy conservation and the hot equation of state 
The energy conservation law gives 



r(r-i) 

2A 



- 1 



p (i - /i)r - i 



(3.8) 



and hence 



/i = l 



r + 1 i / r + i 



4A 



A V 4A 



p 



l + r 



r+i i /r+i 



(3.9) 
(3.10) 



4A y A \ 4A 

where the sign of the square root is the correct one for the Lax shock. In the limit A — > 0, 

T - 1 



/' 



r + i 



and 



p, _ r(r + i) 
p " 



2A 



(3.11) 



(3.12) 
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which is equivalent to 




(3.13) 



This last result means that in the limit in which the internal energy in front of the shock 
can be neglected, the internal energy behind the shock is just the kinetic energy of the 
impact. 



This is the same limit as in the cold case, because for weak shocks shock heating is 
negligible. 

4. The solution behind the shock with the cold EOS 



Tracing a Rieman invariant into the shock from behind, one sees that with c and t finite, 
and the initial conditions finite, v s must remain finite. Generically it will not be zero. As 
the shock approaches the surface of the star, where cq — > 0, we have A — > 0, and hence 
H -> 0. 

From p.4[) we then have v s < a. It is natural to assume that v s has a limit v s — > v* as 
the surface is approached. This parameter v» is determined by the initial data. Therefore, 
near the surface we have v s ~ a ~ «*, and hence, after adjusting the origin of t, 



In the opposite limit A — ¥ oo, which will not be relevant here, we have fi ~ 1 — A 1//2 . 



4.1. Matching conditions at the shock 



(4.1) 



(Note that in our conventions x s < and t < 0.) From (|2.17[) , we have 




gx s 



(4.2) 



n 



n 



and hence 
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From this with ()3.7[) we have 

V*(-zY> ( 4 - 4 ) 



r 



where we define 



Hence we find 



riTv* 
t = 



£ = v*t. (4.5) 



c s = c /i~^ T 1 / 2 ^ (-Vj ^ . (4.6) 
To the same approximation (leading order in /i) , we can write Q3.2p as 

a — v s ~ (j,v* ~ f ^ . (4.7) 

Eqs. (|4.6|) and (|4.7|l provide boundary conditions at the shock for the smooth solu- 
tion behind it. Our derivation of Eq. (14. 6p generalises the "Whitham approximation" of 



H. B. Keller & Whitham(1959) from the shallow water case n = 1 to general n. 

We see that Kq drops out of the result when expressed in terms of c rather than p, 
but both v* and g are relevant. (Note that g comes in not through its effect on the fluid 
behind the shock, which we neglect, but from the assumption that the fluid in front of 
the shock is in hydrostatic equilibrium.) 

4.2. Similarity solution 

The exact solution behind the shock is likely to be attracted to a similarity solution, 
with the initial data remembered only through the parameter v* . To remove v* as far as 
possible, we go to a frame moving with constant velocity , and define the new variables 

u = v — v*, £ = x — v*t (4.8) 

The PDEs that apply behind the shock are Galileo- invariant, and so (|2.11j) becomes 

[d t + (u± c)i%] (u ± 2nc) +.g = 0, (4.9) 
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but in the shock conditions we must take into account that the unperturbed fluid has 



v = 0. Therefore we also define 



^=a-v*, (4.10) 



and note 



s — u s = a — v a . (4-11) 

It is clear that in the limit t — > 0, we can neglect the gravitational force compared to 
pressure forces, and so we look for a solution with g = 0. 
We can write the similarity solution as 

c(x,t) = U(y) (4.12) 



u(x,t) = U(y), (4.13) 



where 



I \ (4-14) 



-K 


t 




T 



with (3 still undetermined. (See Appendix 151 for a derivation.) 

Above, we have shown that the matching conditions at the shock require 

c s ~t*, (4.15) 

s - u s ~ < c s (4.16) 

as f — > 0_, and so 

s~u s . (4.17) 
From the last equation, comparing with (|4.13p . we have y s ^ const to leading order. The 
matching condition (|4.17[) gives 

u{y s ) = p. (4.18) 

Note that the more accurate matching condition (|4.7[) could only be imposed at a higher 



Universality in the run-up of shock waves 15 
order in A as A — > 0. The matching condition (|4.6|) then gives 



y s c{y s ) = V 1 ' 2 , (4.19) 

„ 2 + 3n 2r + 1 . 
"=2 + 2^ = ^ (4 - 20) 

for the overall factor and exponent, respectively. 

The relevant similarity solution, for arbitrary /3, can be found in closed form under the 

assumption that u + 2nc = 0, which means that the Riemann invariant running into the 

shock from behind vanishes (see Appendix [Bj , or physically, that the solution behind the 

shock is dominated by a wave travelling backwards from the shock. This solution can be 

found in implicit form as 



i-0 A 



(4.21) 



y 

u{y) + 2nc(y) = 0. (4.22) 



The matching conditions (|4.18[) and (|4.19|) then fix 



A=Vr( 2 + 3n + 2n lY 2+2 \ (4.23) 
\2 + 7n + 6n 2 J ' V ' 

2 + 3n 



We can write 

-l 



T 



c(x,t) = u» ( -- 1 c(y), c{y) = yc{y), (4.25) 



and similarly for u(x, t). Hence c(x, t) is a regular function of x at t < if and only if c(y) 
is a regular function of y. From (|4.21[) the latter is the case for the range j/ m ; n < y < oo 
for some y min < 0. 

As discussed in Appendix [B] the sonic point y = y c is a line of constant y which is 
also a fluid characteristic. It is given by the condition c = —fi/(2n — 1). Fig. [T] proves 
graphically that y mm < y s < y c < oo for all n > 0. (The expressions for y m j n and y c are 
known in closed form but are long.) Therefore, for t < the similarity solution exists 



16 C. Gundlach and R. J. LeVeque 




Figure 1. Cold EOS: Plots of (from below) j/ m i n (lower boundary of the domain of the similarity 
solution, red, dashed), y s (shock location, blue, dashed) and y c (sonic point, green, solid) against 
the polytropic index n. 

between the shock location x = x s and x = — oo, passing through a sonic point x = x c 
closely behind the shock. Fig. [2] shows c(y) against y for n = 1. 
We can also write 

c(x,t) = (-£) c(y), c{y) = yh(y), (4.26) 



f 



and similarly for u(x,t), where 



^ 1 -^ = 2T^ = 2TTI- (4 - 27) 

The limit 

c(oo) = A 1 " 3 {2n+ l)-~< (4.28) 
exists, and so u(x, t) and c(x,t) are regular as t — > 0_, and become simply 

c(x,0) = v^oo) (-^) 7 , u(x,0) = -2nc(x,0). (4.29) 

4.3. Evolution after the shock has reached the surface 

The data (|4.29[) at t = are smooth, and their evolution to t > remains smooth. We 
were able to neglect gravity in the approach of the shock to the surface, but we must 
restore it in the subsequent smooth motion which occurs on a much longer timescale. 
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Figure 2. Cold EOS: Plot of c(y) against y over the range y a < y < 10, for n = 1. This plot 
gives the instantaneous profile of the sound speed, compare Eqs. (|4.12[) and (I4.14|l . 
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Table 1. Cold EOS: Table of approximate numerical values (the closed form expressions are 
given in the text) of j3, 7, A, y s , y c and c(oo) against the polytropic index V. 



This can be done exactly by going to a freely falling frame by replacing u and £ with 

C = x + vJ-^gt 2 , (4.30) 
w = v — + gt. (4-31) 
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The equations become 

[d t + (w± c)d c ] (w ± 2nc) = 0. (4.32) 

Clearly w + 2nc = at t = and hence for t > 0, and so the solution for t > is a simple 
wave characterised by 

ct - (2ra + l)cc c = 0. (4.33) 
The method of characteristics gives c(£, t) implicitly as 

c = / [C + (2n + l)rf] , (4.34) 

where /(a;) = c(x, 0) given in (|4.29|) . 

5. The solution behind the shock with the hot EOS 

5.1. Matching conditions at the shock 

Our numerical experiments show that the fluid velocity v s behind the shock does not 
remain constant but blows up as the shock approaches the surface. This implies that 
once again we are in the regime A —> 0. From p. 21) with p. Ill) we obtain 

cr~— ^— V s . (5.1) 

With (|2~TU|) . p~T3"l) gives us 

We need a third matching condition for the density. From (|3.11j) . we have 

Ps - fr—[Pa- (5-3) 
Here the equilibrium density in front of the shock po is given by Eq. (|2.19[) . 
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5.2. Similarity solution 

The similarity ansatz for the velocity and sound speed is similar to (|4.12j|4.f4")) . namely 

c(x,t) = ^c(y), (5.4) 

X 

v(x,t) = -u(y), (5.5) 

p{x,t) = C p {-x) n *~p(y), (5.6) 

where 

y = C y (-x)(-t)-^, (5.7) 

where C y is a parameter (of dimension L~ 1 T@) that survives from the initial data into 
the self-similar regime. 

The relevant similarity solution begins at y = oo (corresponding to x — — oo at constant 
t < 0), goes through a regular sonic point y = y Cl and ends at the shock at y = y s . As 
the similarity equations are autonomous in r\ = lny, we can assume that the sonic point 
is at y = 1. 

We expand the solution through the regular sonic point as u — uq + uin + 0(r/ 2 ) and 
c = Co + c\j] + 0(r] 2 ). For given n*, /3 and n, we obtain a unique value of (uq,Cq), but 
two values of (ui, c\). These correspond to two smooth solution curves going through the 
same regular sonic point. The correct solution is easily identified by plotting. 

With the shock at ?y = r/ s , the matching conditions at the shock give 

" (7?s) = fTT (5 ' 8) 

KVs) = ^j- (5.10) 
The requirement that for given n* and n the solution curve (u(rj) , c{rf)) goes smoothly 
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through the sonic point and ends at the shock then determines (3 and rj s < through a 
boundary value problem. 

We solve this boundary value problem numerically by shooting from the sonic point, 
or rather from a small negative value of r\ reached by the power-law expansion around 
77 = 0, to a trial value of rj s , using a trial value of /3. An accurate first guess for rj s and /? 
is obtained by approximating the solution curve from the sonic point to the shock as a 
linear function of r\ using u\ and c\. 

As in the cold case, we can write 

c(x,t) = C- 1 (-tf- 1 c(y) (5.11) 
= Cy?{-xy~c{y), (5.12) 

and similarly for u(x,t). The limits c(oo) and u(oo) exist and are tabulated in Table 01 
In contrast to the cold case, 7 < 0, and so c and v behind the shock diverge asi->0_. 
Similarly, the instantaneous values of c and v just behind the shock diverge because 
[3 — 1 < 0. In physical terms this can be thought of as a whiplash effect, where a finite 
amount of kinetic and internal energy is injected into ever less mass. 
We can write the instantaneous density profile at t < as 

p(x, t) = C p C- n * (-tr>Pp(y), p(y) = y n *p(y). (5.13) 

Finally, the limit p(oo) exists and is tabulated below. Therefore the density immediately 
behind the shock, and the density everywhere as t — > 0_ remain finite, as one would 
expect from the fact that the shock conditions show a finite compression ratio. The 
density profile at t = is just the one in the rest state multiplied by p(oo). This constant 
can be found explicitly in terms of other known constants by evaluating the integral (jC 5[) 
of the similarity equations at y = y s and as y — > 00. 



As a test of our calculation, Table [2] reproduces Table 1 of Sakurai(1960) Focussing 
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r\n* 



1/2 



5/3 
7/5 
6/5 



0.435628 0.22336 0.1141 

0.3934 0.202151 0.103519 
0.330985 0.170418 0.0874924 



Table 2. Nonrelativistic hot EOS: Recreation of Table 1 of Sakurai(1960) 1//3 tabulated 

against F (down) and n* (right). 



n\n* 



3/2 



3 0.642161 
3/2 0.608205 0.751791 
1 0.592366 0.739586 0.807808 
Table 3. Nonrelativistic hot EOS: Table of f3 against n and n*. We only consider stable 

stratifications where n, > n. 



then on the case of an isentropic rest state, Table U gives key parameters of the similarity 
solutions for selected values of T. 

5.3. New hydrostatic equilibrium after the shock 
Consider now an arbitrary instantaneous state in which density and internal energy follow 
approximate power laws, 

p~x a , (5.14) 
cl~T~e~ x b . (5.15) 

Note that at t = in the shock heating process we have discussed, b = 27 and a = n*. 
By contrast, Eq. (|2.21|) tells us that in hydrostatic equilibrium 6=1. 
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p 




ft 

p 










p(oo) 


1.1 





436716 


-0.0211101 





127454 


0.0190103 


765.341 


1.2 





555056 


-0.0468527 





303072 


0.0664771 


2.01223 


1.3 





624385 


-0.0722611 





356576 


0.101789 


1.11067 


4/3 





642089 


-0.0802736 





368456 


0.113245 


0.994484 


1.4 





672222 


-0.0954616 





385789 


0.135492 


0.864565 


1.5 





707991 


-0.116087 





400758 


0.167296 


0.808397 


1.6 





736085 


-0.134239 





410263 


0.198712 


0.743478 


5/3 





751778 


-0.14509 





413301 


0.218764 


0.363108 


1.7 





758898 


-0.150172 





414418 


0.228711 


0.673344 


1.8 





777879 


-0.164173 





414795 


0.257156 


0.655824 


1.9 





793966 


-0.176514 





412426 


0.283948 


0.63714 


2 





807803 


-0.187436 





408344 


0.309271 


0.62677 



Table 4. Nonrelativistic hot EOS: Table of the parameters ft, n 3 and asymptotic values at 
y — oo of u, c and p, against the polytropic index n, for the isentropic rest state case n* = n. 
r > 2 is not considered because the sound speed at ultrarelativistic temperatures is greater than 
the speed of light. F = 1 is unphysical because the speed of sound is zero. We also have difficulty 
reaching r = 1 when numerically solving the boundary value problem for the similarity solution. 



From (|5.14j) and (|5.15|) we have 

K~ x b -". (5.16) 

We define the total mass per area (measured from the surface inwards) as 

,o 

m= I pdx~x a+1 (5.17) 

J X 

and hence 

K~m^. (5.18) 
Let us now assume that the fluid reaches a new hydrostatic equilibrium starting from 




0.2 0.4 0.6 0.8 1.0 



Figure 3. Nonrelativistic hot EOS: Plot of the similarity solution with n = n» = 1 in the 
uc plane, illustrating the dynamical systems ideas discussed in Sec. 15.21 u is right, and c is up. 
The vertical straight line is the flow line and the diagonal straight line are the sonic lines (see 
Apppendix [Cj . The endpoint at u = c = corresponds to y = oo, while the other endpoint 
corresponds to the shock y — y 3 , with the sonic point at y c = 1 by construction. 




Figure 4. Nonrelativistic hot EOS: Plots of (from above) u(y) (blue, dashed), c(y) (green, dot- 
ted) and p/20 (red, solid) for the similarity solution with n = n* = 1 for the range y a < y < e . 
The sonic point is at y c = 1 by definition. These curves give the instantaneous profiles of u, c 
and p against x at finite t < 0, compare Eqs. (|5.7[) . (|5.11[) (and similarly for u) and (|5.13[) . 

the state (|5. 1415. 15|) . This new equilibrium has a density distribution with power a and 
by definition has b = 1 . If heat conduction, shock heating and other entropy-generating 
processes can be neglected in the approach to equilibrium, and if the fluid motion remains 
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plane (spherically) symmetric, then the function K = K(m) is unchanged, because m is 
a Lagrangian coordinate for the plane-symmetric motion. In particular, we must have 

b-a 1 _ a 



a + 1 a+ 1 

which solves to 



(5.19) 



(n + l)a + n(l-6) 

1 + 6n v ; 



6. Numerical experiments 

6.1. Numerical method 
To confirm the validity of the similarity solutions presented here, we have done extensive 
numerical simulation using a high-resolution finite volume method designed to accurately 



capture shock waves. We use the open-source Clawpack software of LeVeque et aZ.(2010) 



The codes used to produce the figures in this section are available at Gundlach & LeVeque(2010) 
along with more figures and animations of the results over time. 

This software requires a "Riemann solver" as the basic building block. For a homoge- 
neous conservation law with no source term, the Riemann solver takes cell averages in 
two neighboring grid cells and determines a set of waves propagating away from the cell 
interface in the solution to the Riemann problem (the conservation law with piecewise 
constant intial data). A simple update of the cell averages based on the distance these 
waves propagate into the cells gives the classic Godunov method, a robust but only first- 
order accurate numerical method. Second order correction terms can be defined in terms 
of these waves and then limiters are applied to these terms in order to avoid non-physical 
oscillations in the solution. This is crucial near discontinuities for problems involving 
shock waves, particularly when near a vacuum state as in the present problem. Com- 



plete details of the algorithms implemented in Clawpack can be found in LeVeque(2002) 
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LeVeque(1997) Several other books also discuss shock-capturing finite volume methods 



based on Riemann solvers, such as Toro(1997) and Trangenstein(2009) The development 
of high-resolution shock-capturing methods has a long history, and many references can 
be found in the books cited above. 

The discussion here will focus on several modifications that have been made to the 
standard software in order to handle the demands of the present problem. We use 
the f-wave formulation of the the wave-propagation algorithm originally proposed in 



Bale et aZ.(2002) in order to obtain a "well-balanced" method that preserves the steady 



state equilibrium solution in the stationary atmosphere ahead of the shock wave, fol- 



lowing the approach in LeVeque(2010) Maintaining this steady state requires that the 
pressure gradient in (|2.2[) balances the source term —gp. In the f-wave formulation, the 
flux difference between adjacent cells is modified by the source terms and the Riemann 
solver applied to this modified flux difference to obtain f-waves that are used in place 
of the classic waves of the Riemann solution. If the source terms are appropriately dis- 
cretized at each interface, the method will be well balanced in the sense that initial data 
that are in equilibrium will lead to zero-strength f-waves in each Riemann solution and 
hence no modification to the solution. Moreover, small perturbations to an equilibrium 
solution will result in small amplitude f-waves. The limiters are applied to these f-waves 
and much more accurate solutions can be obtained for certain quasi-steady problems 
with this approach than with techniques such as fractional step methods (in which one 
alternates between solving the homogeneous conservation law and applying the source 
terms). Discretization of the source term — gp requires defining an appropriate average 
of the densities from the two cells bordering each interface and we use the approach 



suggested in LeVeque(2010) to achieve a well-balanced method for the atmosphere at 
rest. 
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The Riemann solver for the f-wave method must then take the flux difference between 

two cells, as modified by the source term, and split this vector into waves propagating 

to the left and right that are used to update neighboring cell averages. This is typically 

done by splitting the vector into eigenvectors of some approximate Jacobian matrix such 

as the one produced by a Roe average of the neighboring states. However, near vacuum 

state the standard approach can lead to negative pressures or densities and a breakdown 

of the method. We use a variant of the Siliciu relaxation solver described in Sections 



2.4.4-6 of Bouchut(2004) which is related to the HLLE solver. This solver has been 



modified to work with the f-wave formulation in recent work with Murphy(2009) 



To see the self-similar structure that eventually develops, we need to zoom in by many 
orders of magnitude in the vicinity of x = at the times just before the shock reaches the 
origin. For general initial data specified at some time to it is impossible to determine a 
priori the exact final time tf when the shock will reach x = 0, and it varies depending on 
the mesh spacing and the choice of computational parameters such as the limiter used to 
maintain stability. For the computations presented here we have used the monotonized 



centered (MC) limiter (see LeVeque(2002) ), which is generally a good limiter for nonlinear 
systems that gives sharper results than the minmod limiter but greater stability than 
superbee, for example. We have also tried minmod and obtain similar results. Note that 
the parameter t used elsewhere in the paper for the similarity solution corresponds to 
t = t c — tf where t c is the time variable in the computation. We start the computation 
with t c = but reset to t c = at some point as we approach the final time in order to 
retain more significant figures in t. 

We compute on a domain a; m i n < x < which is divided into M finite volume cells. 
Initially x m i n = —12, but we repeatedly cut the domain in half by halving both a; m i n 
and Ax in order to zoom in on the origin. In each time step we estimate the current 
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shock location x s as the location of the largest change in density and we perform this 

rcgridding as soon as x s > 0.1a; m i n , i.e., when the shock is 90% of the way to the 
boundary. Rcgridding is performed by taking the solution on the right half of the domain 
and doubling the resolution, so that M remains the same while a; m ; n and Aa; are halved. 
In each grid cell j = M/2 + 1, . . . , M we estimate a slope for each component of the 
solution by applying a limiter to the one-sided slopes to the left and right, and use this 
to construct a linear function passing through the cell average with this slope (essentially 
identical results are obtained with either the minmod or superbee limiters). We sample 
the resulting linear function at the midpoint of the left and right halves of the cell. These 
two values replace the values stored in cells i and i + 1, where i = 2(j — M/2) — 1. This 
procedure is conservative since the two new values have the same mean value as the 
previous single value. It is second order accurate where the solution is smooth, but it 
does introduce some artifacts at the shock. The new interpolated values do not lie exactly 
on the numerical shock Hugoniot for a propagating shock and hence small oscillations 
generally emerge from the shock at the rcgridding times that propagate into the flow 
behind the shock. These are quite minor and do not hamper our ability to see the 
expected self-similar behavior. 

The time step At is chosen automatically by the Clawpack software, based on a user- 
supplied target Courant number, which we take to be 0.8. The Courant number is com- 
puted in each time step as 

OF£ = ^maxlsfl (6.1) 

Ax i,P 

where s^ is the wave speed of the pth wave for the Riemann problem between cells i — 1 
and i. The time step At for the next step is chosen by setting CFL in (|6.1|) to the 
target Courant number and solving for At. (The actual Courant number in the next step 
may differ from the target since the nonlinear wave speeds are recomputed in each 
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step. If the actual Courant number is larger than the stability limit of 1 then the step 
is re-taken with a smaller At, which at this point can be calculated to hit the target 
Courant number exactly since the waves speeds are now known.) Using this mechanism, 
the time step automatically scales with Ax as the computational domain is subdivided. 
We output results every N out time steps (some fixed number depending on M ) to obtain 
snapshots of the solution as it gets further into the asymptotic regime. The space-time 
numerical grids used resemble a set of nested boxes all sharing the top right corner at 
x = 0, t c = tf. Note that with this strategy the shock never reaches x = and we 
typically halt the computation when ~ 1CU 10 (after about 35 domain halvings). 

From the final frames of the solution we can estimate tf, the time when the shock 
would hit x = 0. As mentioned above, to obtain enough significant figures in tf — t c we 
reset t c to at some point in the computation, typically when x m i n ~ — 10~ 6 . 

For boundary conditions at x = x m i n we use the zero-order extrapolation conditions 
that are often used in Clawpack to avoid non-physical reflections at computational bound- 
aries. The values in the first interior grid cell are simply copied to ghost cells adjacent 



to the cell before each time step as described in Chapter 5 of LeVeque(2002) This is 
obviously unphysical in that it does not represent the overall solution that we have cut 
off. However, this has no influence on the portion of the solution we are interested in. 
After regridding, the shock has been shifted from x w 0.1x m j n to x ~ 0.2a; m i n after x m ; n 
is halved. Using an explicit method with Courant number less than 1, information can 
propagate no more than 1 grid cell per time step. Hence in the time between regriddings 
any effect of boundary conditions at x = a; m i n can contaminate at most 10% of the cells 
before the next regridding takes place, at which point we throw away the left half of the 
cells. So these boundary conditions can never affect the solution near the shock wave. 
The boundary condition imposed at x = is immaterial for this work since we never 



Universality in the run-up of shock waves 



29 



allow the shock to reach the boundary. We need only insure that the equilibrium solution 
is not disturbed adjacent to the boundary. We use a well-balanced method as described 
above that maintains the steady state between regridding times. Unfortunately, this is 
disturbed by regridding since reconstructing and sampling a piecewise linear function 
as described above gives new cell values that are not in exact equilibrium. So in the 
undisturbed region ahead of the shock we reset the solution by evaluting the exact form 
of the known equilibrium solution on the new finer grid each time we regrid. 

6.2. Cold equation of state 

For the cold EOS with n — 1, the equilibrium solution has p(x) = —x, v(x) = 0, 
and c(x) = \/~x. For numerical experiments we perturbed this with strong generic 
displacements in density and momentum centered about different points, namely 



on — 12 < x < 0. The late-time behaviour is independent of these initial data except for 
the fitting of the free parameters v* and t c , thus showing that the similarity solution is 
universal. Figures [5] through [7] show results for this case. 

For this case, M = 4800 grid cells were used and we took 40000 time steps, at which 
point the shock was within distance 1 x 10 -14 of the origin. This required less than 3 
minutes of CPU time on a MacBook Pro using a single 2.2GHz processor. 

A wave travelling right has sharpened into a double N-wave by the time its leading 
edge has reached x = —8.5. By the time the leading edge has reached x — —3.4, the 
leading N-wave has been caught up by and merged with the second N-wave, and the 
wave going left has steepened into an N-wave and completely left the grid. 

The expected similarity solution is obtained by solving Eq. (|4.21[) , multiplied by y, 




(6.2) 




(6.3) 
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Figure 5. Numerical solution for the cold EOS n = 1.0, showing the fluid velocity and sound 
speed after 0, 800, 1600, and 8000 time steps on a grid with M = 4800 grid cells. The generic 
Gaussian initial data breaks up into left-going and right-going waves. The times listed here and 
in later figures are t c — tf as discussed in the text, corresponding to t of the similarity solution. 




Figure 6. Computed shock location for the cold EOS with n = 1, illustrating that the shock 
has constant speed «* as t c — > tf over many decades. Results of more than 21000 time steps are 
plotted (with M = 4800) on a log-log scale. 

numerically for c(y) = yc(y). The parameters v* and tf (the final time, when the shock 
arrives at the surface) must be fit from the computed solution. For late times in the 
computation, the shock appears to propagate at nearly constant velocity and fitting a 
linear function to the observed shock location gives an estimate for tf. The constant 
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Figure 7. Cold equation of state with n = 1. The density (top), velocity (middle), and sound 
speed (bottom) at two different times. The left column shows the solution after 20000 time steps 
and the right column shows the solution after 24000 time steps. In all cases the solid blue curve 
is the computed solution with M = 4800 grid cells and the red dashed curve is the similarity 



solution. For animations, see Gundlach & LeVeque(2010) 



velocity should be v* , but we have found that a more accurate estimate of v* is obtained 
by examining the Ricmann invariant v + 2nc = = const, which is very nearly constant 
behind the shock in the computed solution. 

As x s —> 0, the instantaneous profiles of c and v in the numerical solution behind the 
shock begin to resemble the similarity solution first in shape and in the shock location 
while an overall factor in c and v agrees only later. By the time the shock is at x s ~ 10~ 3 , 
there is agreement within a factor of 2. The agreement is convincing by the time the 
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shock has reached x s ~ 10~ 6 . This is a strong indication that the similarity solution is 
an attractor, at sufficiently small x s and t c — tf. The relevant dimcnsionlcss number here 
is that x s is about 10 -6 times its value when a shock first formed from generic initial 
data. As x s — > as t c — tf — > 0_, the agreement is very good at late times, as illustrated 
in Figure [51 

Figure [7] shows the results after 20000 and 24000 time steps, illustrating close agree- 
ment with the similarity solution. At later times the agreement is not so good though 
convergence tests show that at any fixed time there is increasingly good agreement as 
the grid is refined. 

6.3. Hot equation of state 

Figures 181 through ITU1 show results for three hot EOS cases: n = n* = 1, n = n* = 1.5, 
and n = 1, n* = 1.5. The initial behaviour of the shocks is qualitatively similar to the 
cold EOS case. The late-time behaviour approaches the similarity solution and is again 
independent of the initial data except for the fitting of the free parameters C y and t c . 
Note the very different scaling of the axes in each case. 

The hot EOS case is easier to compute than cold EOS case and all of these results 
are shown on grids with M = 1200 cells. Finer grids give even closer agreement. With 
M = 4800 the computed solution is indistinguishable from the similarity solutions to 
plotting accuracy. Moreover, the solutions continue to match well out to x s ~ 10~ 15 , 
after starting to agree with the similarity solution very well at t ~ 10 -3 , when x s ~ 0.01. 

7. Conclusions 

We have constructed similarity solutions that describe a shock approaching the surface 
of a star from the interior, in the limit where the shock is very close to the surface and 
can be approximated as planar. We have shown by numerical simulation in the planar 



Universality in the run-up of shock waves 



33 



density at time t = -2.641e-08 



computed soln, xshock = -1.678e-06 
similarity soln, xshock = -1.674e-06 




0.000010 -0.000008 -0.000006 -0.000004 -0.000002 0.000000 



velocity at time t = -2.641e-08 



computed soln, xshock = -1.678e-06 
similarity soln, xshock = -1.674e-06 



-0.000010 -0.000008 -0.000006 -0.000004 -0.000002 0.000000 



sound speed at time t = -2.641e-08 




computed soln, xshock = -1.678e-06 
similarity soln, xshock = -1.674e-06 



density at time t = -9.831e-15 



computed soln, xshock = -1.077e-ll 
similarity soln, xshock = -1.072e-ll 




600 
500 
400 
300 
200 
100 



velocity at time t = -9.831e-15 



computed soln, xshock = -1.077e-ll 
similarity soln, xshock = -1.072e-ll 



sound speed at time t = -9.831e-15 



400 
300 
200 
100 



— computed soln, xshock = -1.077e-ll 
-- similarity soln, xshock = -1.072e-ll 



-0.000010 -0.000008 -0.000006 -0.000004 -0.000002 0.000000 



-7 -6 -5 -4 



Figure 8. Hot equation of state with n = n* = 1. The density (top), velocity (middle), and 
sound speed (bottom) at two different times. The left column shows the solution after 6000 
time steps and the right column shows the solution after 10000 time steps. In all cases the solid 
blue curve is the computed solution with M = 1200 grid cells and the red dashed curve is the 



similarity solution. For animations, see Gundlach & LeVeque(2010) 



approximation that generic initial data approach the similarity solution as a universal 
endstate as the shock approaches the surface. 

The behaviour is very different for cold and hot equations of state. In the cold EOS, 
the fluid behind the shock moves to leading order as a single body, with the matter still 
in front of the shock getting stuck on like insects to the front of a car, and dynamically 
insignificant. The density directly behind the shock becomes ever larger than that in 
front, and the velocity directly behind the shock approaches a constant. 
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Figure 9. Hot equation of state with n — n» = 1.5. The density (top), velocity (middle), and 
sound speed (bottom) at two different times. The left column shows the solution after 6000 
time steps and the right column shows the solution after 10000 time steps. In all cases the solid 
blue curve is the computed solution with M = 1200 grid cells and the red dashed curve is the 



similarity solution. For animations, see Gundlach & LeVeque(2010) 



With the hot equation of state, shock heating becomes essential, and this indicates 
that the cold EOS approximation is completely wrong for shocks near the surface. The 
ratio of densities in front and behind the shock approaches a constant. The fluid velocity 
and sound speed behind the shock diverge, which means that eventually the motion has 
to be treated relativistically. There would then be a (non-self-similar) transition to a new, 
ultra-relativistic similarity solution, where one also has to assume that P = pe/3. This 



ultra-relativistic solution has been investigated by Nakayama & Shigeyama(2005) but 
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velocity at time t = -2.644e-09 



computed soln, xshock = -8.631e-07 
similarity soln, xshock = -8.580e-07 
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Figure 10. Hot equation of state with n = 1.0, n* = 1.5. The density (top), velocity (middle), 
and sound speed (bottom) at two different times. The left column shows the solution after 6000 
time steps and the right column shows the solution after 10000 time steps. In all cases the solid 
blue curve is the computed solution with M = 1200 grid cells and the red dashed curve is the 



similarity solution. For animations, see Gundlach & LeVeque(2010) 



it would depend on the physical circumstances if it becomes relevant before the perfect 
fluid approximation breaks down (and has to be replaced, for example, by kinetic theory) 
at the surface of the star. 

We hope that our solutions provide at least rough approximations to shock formation 
near a stellar surface in some circumstances. They are also good testbeds for codes 
capable of dealing with stellar surfaces. Our simulations are based on the Clawpack 
software but required some modification to standard methods in order to compute into 
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the similarity regime. In particular, we have used a well-balanced approach that maintains 
the hydrostatic equilibrium ahead of the shock to machine precision and a remeshing 
procedure to obtain well resolved solutions down to x s ~ 10~ 15 . 

In an appendix, we have given a self-contained derivation of the similarity solution. 
We stress how similarity solutions for the cold equation of state can be identified with 
isentropic similarity solutions for the hot equation of state, and under what circumstances 
any conservation law in the original fluid equations gives rise to an integral of the motion 
in the ordinary differential equations resulting from a similarity ansatz. 

CG was supported in part by STFC grant PP/E001025/1. RJL was supported in part 
by NSF Grants DMS-0609661 and DMS-0914942, and the Founders Term Professorship 
in Applied Mathematics at the University of Washington. 

Appendix A. General notes on similarity solutions 

A.l. Transport laws 
Consider a transport law of the form 



where v is a velocity, and the dimension of tp is undetermined. Without loss of generality, 
we can write any power-law similarity ansatz as 



tpt + vip x = 0, 



(Al) 



(A 2) 



v = x 



v(y), 



(A3) 



where 



y = xt 



(A 4) 
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(Any more general ansatz is related to this by redefining y, v and 0.) In order to obtain 

an ODE in y alone, we require 

7 =l-~^/?=— L. (A5) 
P 1-7 

The resulting ODE can be simplified by defining 

v(y) = V~' v{y), (A 6) 

and becomes 

{v - fi) y(p' + av0 = 0. (A 7) 

This redefinition is equivalent to defining 

« = -%). (A 8) 

This no longer contains 7, and makes v dimensionlcss. 

If, in particular, a = 7, and we define (p{y) = y 1 ^ip(y), the transport law becomes 

(v - P) yip' + (v - 1)0 = 0. (A 9) 

A. 2. Conservation laws and integrals 
Consider the conservation law 

ipt + F x = 0. (A 10) 

Without loss of generality, the most general power law similarity ansatz that results in 
an ODE is 

ip = x a 0(y), (All) 
F = x a+1 r 1 F{y) 1 (A 12) 

and this ODE is 

- p0' + F' + ^-tlp = 0. (A 13) 

y 
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If and only if a = — 1, this admits an integral of the form 

F — /3tp = const. (A 14) 

In particular, if F = v<p 7 so that we have 

<p t + {v<p) x = 0, (A 15) 

and making the consistent ansatz (|A 8p . we have 

(v-j3)tp = const. (A 16) 

Usually, however, we cannot arrange for a = — 1. 

By contrast, if we have a conserved mass density p, then any quantity K that is 
constant along fluid world lines gives rise to an integral of the ODEs in the self-similarity 
ansatz. From 

K t +vK x = 0, (A 17) 

Pt + (vp) x = 0, (A 18) 

we get 

(pK*) t + (v P K*) x = 0, (A 19) 

P = x a ~ P {y), (A 20) 

K = x"K{y), (A 21) 

and (|A 8|) . we can now choose p = — (a + and obtain 

(v - P) pK p = const. (A 22) 

without any restriction on v. 



for any p. Assuming 
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Appendix B. Cold polytropic gas dynamics 



:-!!) 



B.l. Similarity ansatz 



In gas dynamics u and c both scale as velocities. Therefore, with the ansatz 



(Bl) 



(B2) 



the equations (|4.9[) with g = 0, using (|A 9|) . become 



(u ± c - 1) (u ± 2nc) + (u±c- P)y(u± 2nc)' = 0. 



(B3) 



These equations become automous in the independent variable r\ = lny, thus defining a 
dynamical system in the uc-plane. The direction field is reflection-symmetric about the 
line c = 0, and trajectories cannot cross this line, so that we can assume c > 0. 

We note in passing that a self-similar solution can be obtained also for g ^ if we set 
ft = 2. The dimensionlcss self-similarity variable is then y = x/(gt 2 ), from dimensional 
analysis (self-similarity of the first kind) . 



The general solution for /3 = 1 is the trivial one u = uq, c = cq. 

For (3 ^= 1, the lines c = ±(u — f3) define sonic points, where the lines of constant y are 
tangential to one of the two fluid characteristics in the ict-plane. The abstract dynamical 
system is regular there, with dc/du = =Fl/(2n), but the ODEs in y or 77 arc in general 
singular there, with de/drj ~ 1/e, where e is the distance to the line, and so solution 
curves end there with e ~ ^/j/o — r\. 

However, on each line of sonic points, there is one regular sonic point given by c = 
±w/(2n), where de/dr/ is finite (because both its numerator and denominator vanish). 



B.2. Dynamical system analysis 
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They are 

6=±-?-, (B4) 
2n-V 2n - 1 v ; 

Precisely one of these is in the physical region c > 0, depending on the sign of /3 and 
2n-l. 

There are two smooth solutions through each of these points, with slopes 

dc _ 1 (/?-l)(4n 2 + l) + 4n 

+ 2n 4n [1 + 2n(/3 - 1)] ' 1 J 

A further power-law expansion shows that for the first of these solutions c = ±2nu to 

all orders. One of the two ODEs (|B 3p is then identically satisfied, and the remaining one 

becomes separable when written in c alone and is independent of the sign ± above. An 

implicit solution in closed form is given in (|4.2H4.2"^]l above. 



Appendix C. Hot polytropic gas dynamics 

C.l. Similarity ansatz 
A useful form of the general power-law similarity ansatz is 

P = x n '~p{y), (CI) 
u=ju(v), (C2) 
c=jc(y). (C3) 

In the variable r\ = lay the the resulting ODEs are autonomous. Furthermore, their 
right-hand sides are rational functions of only ii, c and the parameters n», j3 and n. This 
means that we have a dynamical system in the tic-plane, with lnp slaved to u and c. 
Again we assume that c > 0. 
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C.2. Integral of the motion 

The reason that p does not appear on the right-hand sides is that K is constant on fluid 
worldlines, and hence that an integral of the motion exists. We have 

TK = c 2 p-^ = x 2 ~^t- 2 c 2 p-i 

= x 2 ~^~iyic 2 p-" (C4) 

From our discussion in Sec. IA.2I we now have the integral 

(u - P) yfc^p 1 -™ = const., (C 5) 

where 

P= , 9 = 2 ~. (C6) 

q n p 

C.3. Isentropic limit 

To compare the ODEs for u and c with their counterparts for the cold EOS, we can write 
them in the form 

{u ± c - 1) (u ± 2nd) + (u±c~ P)y{u± 2nc)' 

= -^MV (C7) 
1 u — p 

The right-hand side in this equation corresponds to the right-hand side in Eq. (|2.13l) . 

The cold similarity equations (|2.13[) are obtained as the case q = of the hot EOS 
similarity equations (|C 7[) . Conversely, the limit q —> of the integral (jC 5p is 

y%c(y) 2 p(y)-K = const. (C 8) 

and from (jC 4j) we see that K is constant. Therefore, any combination of n*, j3 and n 
such that q = represents an instant of the isentropic limit of the hot EOS case, and 
conversely, any similarity solution with the cold EOS can be interpreted as an isentropic 
similarity solution with the hot EOS, with n* given by q = 0, and p(y) given by (|C 8[) . 
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C.4. Dynamical system analysis 

In general, u'(r)) and c'(j]) are singular (their denominators vanish) on the straight lines 
u — (3 = and c = ±(u — (3). These lines divide the half-plane c > into four sectors. 
Limits y — > and y — > oo 

For |c|, |u| S> 1, \(3\, we have u,c ~ as y — >• and and for |c|, \u\ <?C 1, \/3\, we have 
u, c ~ y~ x ^ as y — > oo. The limit of m/c is undetermined. 
Flow line u = j3 

The line u = (3 represents points where a fluid world line coincides with a similarity 
line. Let us call it the flow line. Near it, u — ~ 7? — 7?o an d c ~ (ry — ^o) ' 39 ^ 21 "- So 
solutions end there at finite ?7 but in the dynamical system it is an asymptote. 

Sonic lines c = ±(u — 0) 

The two other lines represent sonic points, where a sound characteristic coincides with 
a similarity line of constant y. We shall call them the sonic lines of the dynamical system. 
The dynamical system is regular there, with dc/du — =pl/(2n), but the ODEs in y or rj 
are singular there, with de/drj ~ 1/e, where e is the distance to the line, and so solution 
curves end there with e ~ \/?7o ~ V- O n each sonic line, there is also one regular sonic 
point where the y or rj derivatives are also finite. Exactly one of these is in the physical 
region c > 0, depending onti„ /3 and n. Through each regular sonic point there are two 
smooth solutions. 

C.5. Special case (3 = 1 

For (3 = 1, the regular sonic points move down to c = 0, u = 1. There the solution ends 
because c cannot change sign. 

For (3 = 1 only, the Galileo transformation u — >• it + relates different solutions of 
the similarity ODEs, so that any solution (u(jj),v(jj)) can be obtained from a reference 
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solution (u*(y), v*(y)) as 

.4 
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u(y)-l = ( 1- - ) &.(y-A)-l 



(C9) 
(CIO) 

where A is a real constant. Note that this transformation does not mix solutions in the 
four sectors, so one needs one reference solution in each sector to fill the phase plane. 
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